#include "Tria6NodeShpFunc.h"
#include "DataStructure/Common/Real.h"

using SG::DataStructure::Common::Real;
using SG::Algebra::Matrixd;
using SG::FEMSolver::Element::Tria6NodeShpFunc;

Tria6NodeShpFunc::Tria6NodeShpFunc ()
:NSet (), DNDxiSet ()
{}

/** @brief 计算积分点处形函数及形函数对参数坐标导数
 * 
 */
void Tria6NodeShpFunc::Compute (const SG::DataStructure::FEM::IntegPointArray& points)
{
    const std::size_t nGauss = points.size();
    NSet.resize(nGauss);
    DNDxiSet.resize(nGauss);

    for (std::size_t i = 0; i < nGauss; ++i)
    {
        const auto& xi  = points[i].m_x;
        const auto& eta = points[i].m_y;
        const auto  L1  = 1.0 - xi - eta;
        const auto L2 = xi;
        const auto L3 = eta;

        // ---- 形函数 ----
        NSet[i]       = Matrixd (1, 6);
        NSet[i](0, 0) = L1 * (2.0 * L1 - 1.0);
        NSet[i](0, 1) = L2 * (2.0 * L2 - 1.0);
        NSet[i](0, 2) = L3 * (2.0 * L3 - 1.0);
        NSet[i](0, 3) = 4.0 * L1 * L2;
        NSet[i](0, 4) = 4.0 * L2 * L3;
        NSet[i](0, 5) = 4.0 * L3 * L1;

        // ---- 形函数对 ξ、η 的导数 ----
        DNDxiSet[i] = Matrixd (2, 6);

        // 对 ξ
        DNDxiSet[i](0, 0) = -4.0 * L1 + 1.0;
        DNDxiSet[i](0, 1) = 4.0 * L2 - 1.0;
        DNDxiSet[i](0, 2) = 0.0;
        DNDxiSet[i](0, 3) = 4.0 * (L1 - L2);
        DNDxiSet[i](0, 4) = 4.0 * L3;
        DNDxiSet[i](0, 5) = -4.0 * L3;

        // 对 η
        DNDxiSet[i](1, 0) = -4.0 * L1 + 1.0;
        DNDxiSet[i](1, 1) = 0.0;
        DNDxiSet[i](1, 2) = 4.0 * L3 - 1.0;
        DNDxiSet[i](1, 3) = -4.0 * L2;
        DNDxiSet[i](1, 4) = 4.0 * L2;
        DNDxiSet[i](1, 5) = 4.0 * (L1 - L3);
    }
}